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PROGRAM SUMMARY 

Title of program: SABSPV 

Computer: DEC VAX, HP/APOLLO, IBM/RS6000; Installation: INFN, Sezione di 
Pavia, via A. Bassi 6 , 27100 Pavia, Italy 

Operating system: VMS, UNIX 

Programming language used: FORTRAN 77 

Memory required to execute with typical data: 40 kByte 

No. of bits in a word: 32 

The code has not been vectorized 

Subprograms used: random number generator RANLUX [1]; routines from the CERN Pro¬ 
gram Library are also used. 

No. of lines in distributed program, including test data etc.: 1700. 

Correspondence to: 

CACCIARI@DESY.DE, CACCIARI@PV.INFN.IT; 

MONTAGNA@PV.INFN.IT; 

NICROSINI@VXCERN.CERN.CH, NICROSINI@PV.INFN.IT; 
PICCININI@PV.INFN.IT 

Keywords: high energy electron-positron collisions, small-angle Bhabha scattering, lumi¬ 
nosity determination, QED corrections, electron structure functions, experimental cuts, 
Monte Carlo integration. 

Nature of physical problem 

The precise determination of the theoretical Bhabha scattering cross section in the small- 
angle regime is a key ingredient for precision luminometry at LEP. To this aim, QED 
radiative corrections to the tree-level Standard Model Cross section have to be taken into 
account, together with vacuum polarization effects. Particular care has to be devoted to 
higher-order corrections. The theoretical formulation must allow the computation of the 
cross section for a wide variety of cuts on the final-state particles. 

Method of solution 

A suitable matching of a fixed-order calculation with the structure-function techniques for 
resumming large initial- and final-state leading logarithmic corrections [ 2 ] is performed. 
A Monte Carlo integration with weighted events has been implemented in order to mimic 
as close as possible the experimental triggering conditions. The importance-sampling 
technique [3] is employed to take care of the peaking behaviour of the integrand. 

Restrictions on the complexity of the problem 

The O(a) QED corrections are computed exactly only for the dominant contribution to 
the small-angle cross section, namely the square of 7 -exchange in the t-channel; all the 
other contributions to the cross sections are corrected at the leading-logarithmic level. 
The contribution of additional hadronic or leptonic pairs is at present neglected. Starting 
from 0(a 2 ), QED corrections are implemented at the leading-logarithmic level. 
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Typical running time 

On a HP 9000/735 the code takes about 2 x 10 -4 seconds per event, with standard cuts. 
A one per mille relative error can be achieved with about 10' events. 

Unusual features of the program: none 
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LONG WRITE - UP 


1 Introduction 


Precision luminometry at LEP requires the best possible knowledge of the theoretical 
Bhabha scattering cross section in the small-angle regime. Actually, the experimental 
precision of luminosity measurements is at present at the level of 0.1% and is foreseen to 
improve further in the near future. In order to exploit such an experimental precision, 
theoretical results of comparable accuracy are required. 

Two kinds of problems are to be faced here. The first one is to establish the physical 
precision of a given theoretical approach; this requires an estimate of those contributions to 
the cross section which are neglected by the approach itself. The second one is to establish 
the technical precision of the computational realization of a given theoretical approach; 
this requires an estimate of the errors introduced by numerical approximations, different 
implementations of theoretically known effects and so on, and, last but not least, bugs in 
the computer code. 

The comparison between different and independent approaches can contribute to the 
solution of both problems and become important in order, on the one hand, to assess the 
presently reached accuracy and, on the other, to eventually improve it in a significant 
way. 

Some codes that compute the small-angle Bhabha scattering cross section already 
exist II- The aim of the present work is to describe in some detail the code SABSPV . 

In order to compute a small-angle Bhabha scattering cross section with a precision of 
the order of 0.1% and within realistic experimental cuts, the following tasks have to be 
accomplished: 

1 . to use a theoretical formulation capable of taking into account the effects of multiple 
soft and/or collinear photon emission from initial- and final-state particles, but also 
of providing an exact evaluation of at least one hard photon emission; 

2 . to perform an implementation of the integration within the phase-space limits such 
that a wide variety of experimental triggering conditions can be closely reproduced. 

The first point is solved by developing a formulation, described in which allows the 
matching of a fixed-order perturbative calculation with the structure-function technique. 
The outcome of this matching is a calculation which correctly describes one hard photon 
emission but also allows for the evaluation of multiphoton emission in the initial and final 
state. 

As for the second point, the cuts are implemented by using a Monte Carlo integra¬ 
tion. The four-momenta of the outgoing particles (electron, positron and photon) are 
constructed and fed to a cutting routine that reproduces the experimental triggering con¬ 
ditions. 

In the following we will first review the theoretical background and its implementation 
in the code. The structure of the integrand demands the importance-sampling technique 
to be applied in order to have an efficient Monte Carlo integration. We will discuss the 
way this has been done. Finally, the program structure will be briefly presented and the 
required inputs and the corresponding output will be described. 
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2 Theoretical Background 


The structure-function approach allows the leading-log-corrected cross section in the lab¬ 
oratory frame to be written, taking into account all-orders photonic radiation, in the 
following form: 

4 l } = / dndx 1 dx 2 D {oo \x ll Q 2 )D (oo \x 2 1 Q 2 )J(x 1 ,x 2 ,'d) 

x x 2),i{ x h x 2, Q 2 )i 7 '^ oo) (^2) Q 2 )e(cuts). ( 1 ) 

The detailed description of Eq. (|lj) can be found in Q, f|]. For convenience, we recall that cr 
is the Born cross section, the D's are the structure functions for initial-state radiation and 
the F’s describe final-state photonic emission; ©(cuts) represents the rejection algorithm 
for implementing the experimental cuts. 

It was prompted in the Introduction that when aiming at a per rnille accuracy we 
have to rely also on a complete 0(a) calculation of small-angle Bhabha scattering cross 
section. The complete squared matrix element of the (gauge-invariant) f-channel diagrams 
associated with (hard) photon radiation from the electron line has been computed in j5|. 
Its full expression, including all 0(m 2 ), 0(m 4 ) and 0(m 6 ) terms, reads as follows: 

M 0 + m 2 M 2 + m 4 M 4 + m 6 M 6 | , (2) 

where the single terms are explicitly reported in 0. From this matrix element the 
0(a ) cross section can be computed via Monte Carlo integration as 

of (fc 0 ) = fd[PS\ F(JM _| 2 + |M + | 2 )@(cuts) , (3) 

where (/[PS 1 ] is the phase-space volume element, |M _|-| 2 are the squared amplitudes for 
the electronic and positronic radiation, F is the proper flux factor, and k 0 is a minimum 
energy fraction of the radiated photon, defined as ko = EF n /E. 

The matching between the all-orders leading-log cross section of eq. ® and the O(a) 
one (|3|) takes place in the following way: the order-a content of the leading-log f-channel 
y-exchange cross section is extracted by employing the O(a) expansions of the structure 
functions entering the master formula (]T|), thereby yielding crff,. Denoting by cf +v (k 0 ) 
the f-channel 7 -exchange cross section including virtual corrections plus soft photons of 
energy up to F 7 = k^E (see for instance |J), the cross section in the luminometry region 
can finally be written as 

O-A = - 0'S, 7 + + O-f (fc 0 ) , (4) 

where, we recall, crjf'* is th e all-orders leading-log-corrected full Bhabha cross section 
of Eq. o"ll, 7 is U P t° @( a ) leading-log-corrected cross section, limited to the 

t-channel 7 -exchange contribution, and cr^(ko) is the radiative Bhabha cross section of 
Eq. (§. 

The difference crjf ^ ~ 7 contains the Born contribution of all the Bhabha channels 

but the f-channel 7 -exchange one, plus their leading-log QED corrections resummed to 
all orders, plus the higher-order leading-log QED corrections to the f-channel 7 -exchange 


\M\ 2 = 


2(47tq:) 
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contribution starting from 0(a 2 )] the exact up to 0(a) contribution for the f-channel 
7 -exchange term is then supplied by a^ +v (k 0 ) + atf (k 0 ). 

Equation (jlj) is in the additive form. A factorized form can also be supplied. It has 
the same 0(a ) content but also leads to the so-called classical limit, according to which 
the cross section must vanish in the absence of photonic radiation. It reads 

ctf — (1 + C'vlVlzA 


r H 

1 Y / — 


a; 


7 +V ’( fc o) + <r~(ko) 


o 


(«) A a ) 

LL, 7 J NL ,7 


cr 


a 


(5) 


cr being the Born cross section; C$ L contains the non-log part of the O(a) ) cross 

section, represented by crffi L . A detailed discussion on the physical content of eqs. ([§) 
and ([5]) will be given elsewhere. 


3 Importance Sampling 


Both the leading-log cross section and the real hard photon one present singularities which 
have to be properly treated with importance-sampling techniques to ensure a reasonably 
fast convergence of the Monte Carlo integration. This means dividing the integrand by a 
so-called weight function and redefining the integration variable: 



with y = f p(x)dx. The weight function p(x ) must be chosen such that the ratio f(x)/p(x) 
is smooth enough to produce a small variance. 

Let us start with the leading-log cross section. The structure function D i '°°\x 1 Q 2 ) 
reads 


D ( '°°\x, Q 2 ) = 


exp {i/>(!-7*)}.l ^.1 


r(i +1/3) 

+—/3 2 — 4(1 + x) ln(l — x) + 3(1 + x) lux 
32 


-4 


In x 
1 — X 


— 5 — x 


7T 


ln (S, 


- 1 


(7) 


and displays a singularity in x — 1. This is disposed of by employing, in connection with 
each structure function in the integrand, a weight function that resembles the singular 
part of the structure function itself: 

P(x) = ^P(1~ x)^- 1 . ( 8 ) 

The relation between the original integration variable x and the flat random variable 
r G ( 0 , 1 ) is then 

x = 1 — (1 — r) 2 ^. (9) 

The real hard photon cross section evaluation must instead face the “singular” (in 
fact, only extremely peaked) distribution of the relative angle 0 7 between the photon 
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and the emitting particle when approaching the forward collinear region. Four such 
“singularities” are present, related to initial- and final-state emission from electron and 
positron respectively. The initial-state ones are smoothed by choosing the weight function 

p(cos d^) — (l — - ml/E 2 cos d^j , (10) 

m e and E being the electron mass and energy in the lab frame, respectively. The final- 
state ones are disposed of by using similar weight functions after a proper rotation of the 
reference frame such that the z-axis is taken along the outgoing lepton. 


4 Program Structure 


The code has the usual structure of Monte Carlo programs, with routines evaluating the 
cross sections (some pertaining to the real hard photon part, some to the leading-log 
cross section), triggering routines, determining whether the event is to be accepted or 
not, utility routines, and finally routines collecting the output. More in detail, we have: 

• Monte Carlo routines: 

- the MAIN PROGRAM: where the Monte Carlo loop is located. It provides the ran¬ 
dom numbers to the routine INTEGR, collects its outputs, performs the statistics 
operations and writes the results to the output hie; 

- SUBROUTINE INTEGR: receives the random numbers from the MAIN PROGRAM 
and calls the cross sections, kinematics and triggering routines. It returns the 
cross sections for a given event to the MAIN PROGRAM. 

• Cross section routines: 

— SUBROUTINE HARD, SUBROUTINE EEG: provide the real hard photon cross sec¬ 
tion; 

- FUNCTION BHABORN, FUNCTION BHAB0RN_gtgt, FUNCTION AINTSV: provide the 
full Born cross section, the Born cross section for the 7 (f)y(t) channel and the 
virtual plus soft 0(a) correction to the latter, respectively; 

- FUNCTION DOP, FUNCTION AJAC, FUNCTION PDINFTY, FUNCTION PDPDALPHA: 

routines related to the structure-function evaluation of the leading-log cross 
section; 

- FUNCTION VPOL, FUNCTION SOMMAT, FUNCTION ABC: evaluate the vacuum po¬ 
larization correction to the photon propagator, according to the parametriza- 
tion of [ 0 . 

• Kinematic routines: 

- SUBROUTINE KINE2, SUBROUTINE QUADRI2: build up the outgoing particles 
four-momenta from the generated random numbers for two-body events, i.e. 
Born and virtual plus soft cross sections, and collinear leading-log ones; 

- SUBROUTINE KINEH, SUBROUTINE QUADRIH: the same but for three-body events, 
with the real hard photon. 
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• Triggering routines: 

— SUBROUTINE CUTS: master cut routine. It propagates to the rest of the code 
the result of one of the triggering routines; 

— SUBROUTINE TRIGGER1: performs acceptance cuts for a sample trigger; 

— SUBROUTINE TRIGGER2: dummy routine to be used to implement a new trigger. 

• Utility routines: 

- FUNCTION EXCHANGE, FUNCTION COSANGLE, FUNCTION PHIROT: they exchange 
any two variables, evaluate the cosine of the angle between any two four- 
momenta and rotate a four-momentum around its ^-axis, respectively. 

The code can be used in a “black box” way with the only exception of the triggering 
routine, which may be provided by the user. An example of how to build such a routine 
is given in Section |5]. We now want to describe the common blocks the user has to be 
familiar with for writing his own triggering routine. 

• REAL*8 QUADRIM(0:3), QUADRIP(0:3), KUADRI(0:3) 
COMMON/MOMENTA/QUADRIM,QUADRIP,KUADRI 

This common block contains the four-momenta of the outgoing electron (QUADRIM), 
positron (QUADRIP) and photon (KUADRI). It is filled by the QUADRI2 or QUADRIH 
routines when treating soft-collinear or real hard photons respectively. 

• COMMON/CUTOFF/EGMIN,IHARD 

These are non-physical parameters. EGMIN contains the value of the cutoff k o, which 
separates virtual plus soft and real hard photon events. IHARD is a flag which is set 
to 1 when processing real hard photon events, and 0 otherwise. It can be used in 
the triggering routines to skip or implement some cuts according to the event being 
of two-body/collinear type or with a real hard photon. For example, the cut on the 
ko cutoff only has to be activated in the real hard photon part, being taken into 
account analytically in the virtual plus soft part. 

• C0MM0N/C0NST/EBEAM,ALPHA,CONVFAC,PI,AME 

Some constants, which are used in many points of the code. EBEAM is the beam 
energy in GeV (see Section about input parameters). ALPHA is the electromagnetic 
coupling, a = 1/137.0359895. CONVFAC is the conversion factor from GeV -2 to nb: 
.38937966 x 10 6 . Pi is 7r and, finally, AME is the electron mass: m e = 0.51099906 x 
10- 3 GeV. 

• C0MM0N/EXPCUTS/T1MIN, T1MAX, T2MIN, T2MAX, E1MIN, E2MIN 

The experimental acceptance cuts. See Section | about inputs description. 

• REAL*8 CAL0INPUT(5) 

C0MM0N/CAL0S/CAL0INPUT 

These are parameters which are read in from the input hie (once more, see Section |6|) 
and which can be used in writing the triggering routines. The meaning of the five 
parameters will depend on the use which is made of them in the routines. For 
instance, in TRIGGER1 the parameter CALOINPUT(l) is the calorimetric threshold as 
a fraction of the beam energy and CAL0INPUT(2) is the half-opening (in radians) of 
the cone, which defines an electromagnetic cluster. 
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5 Sample Trigger Routine 


A sample triggering routine TRIGGERlO is inserted in the code. It can be used as an 
example to construct other experimental-like triggers. A dummy TRIGGER2 0 is also 
given, to be completed by the user. 

Below we give the listing of TRIGGERlO. It first checks for the electron and the 
positron to be within angular acceptance cuts. After that, it checks whether one of the 
leptons and the photon do form an electromagnetic cluster, i.e. are radially separated by 
less than the critical distance fixed by the input parameter DEL. If they do, the energy cut 
is performed using the sum of their energies rather than the energy of the bare fermion 
only. 

LOGICAL FUNCTION TRIGGERlO 
IMPLICIT REAL*8 (A-H.O-Z) 

REAL*8 QUADRIM(0:3), QUADRIP(0:3), KUADRI(0:3) 

REAL*8 ZETA(0:3) 

REAL*8 CAL0INPUT(5) 

COMMON/MOMENTA/QUADRIM,QUADRIP,KUADRI 

COMMON/EXP CUTS/T1MIN, T1MAX, T2MIN, T2MAX, E1MIN, E2MIN 

C0MM0N/CUT0FF/EGMIN,IHARD 

C0MM0N/C0NST/EBEAM,ALPHA,CONVFAC,PI,AME 

C0MM0N/CAL0S/CAL0INPUT 

DATA ZETA/0D0,0D0,0D0,1D0/ 

* 

EGAMMA = KUADRI(0) 

* 

IF(IHARD.EQ.1) THEN 

* 

* .THIS IS THE NON-PHYSICAL CUTOFF 

* 

IF(EGAMMA.LT.EGMIN) THEN 
TRIGGER1 = .FALSE. 

RETURN 
END IF 
END IF 

* 

THETA1 = ACOS(COSANGLE(QUADRIM,ZETA)) 

THETA2 = ACOS(COSANGLE(QUADRIP,ZETA)) 

THETA20PP = PI-THETA2 
El = QUADRIM(O) 

E2 = QUADRIP(0) 

* 

IF(THETA1.LT.T1MIN) THEN 
TRIGGER1 = .FALSE. 

RETURN 
END IF 

* 

IF(THETA1.GT.T1MAX) THEN 
TRIGGERl = .FALSE. 




RETURN 
END IF 

IF(THETA20PP.LT.T2MIN) THEN 
TRIGGER1 = .FALSE. 

RETURN 
END IF 

IF(THETA20PP.GT.T2MAX) THEN 
TRIGGER1 = .FALSE. 

RETURN 
END IF 

EMCL = El 
EPCL = E2 

IF (IHARD.EQ.l) THEN 
DJET = CAL0INPUT(2) 

CDJET = COS(DJET) 

CEG = COSANGLE(QUADRIM, KUADRI) 
CPG = COSANGLE(QUADRIP, KUADRI) 

IF(CEG.GE.CDJET) THEN 
EMCL = El + EGAMMA 
END IF 

IF(CPG.GE.CDJET) THEN 
EPCL = E2 + EGAMMA 
END IF 

END IF 

CALOTH = CALOINPUT(l) 

SUMEN = EMCL*EPCL/EBEAM/EBEAM 

IF(SUMEN.LT.CALOTH) THEN 
TRIGGER1 = .FALSE. 

RETURN 
END IF 

TRIGGER1 = .TRUE. 

END 


6 Input Description 


A data card of the following kind has to be provided, via standard input, when running 
the program: 
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46.15D0 

24.D-3 58.D-3 O.DO 

24.D-3 58.D-3 O.DO 

0.5D0 l.D-2 O.DO O.DO O.DO 

1 
1 

1.D5 O.DO 0 ’ SABSPV.OUT’ 


EBEAM 

T1MIN, UMAX, E1MIN 
T2MIN, T2MAX, E2MIN 
CALOINPUT(l...5) 

ISIM 

ICALO 

EVTS, ACCLIM, IRESTART, OUTFILE 


These parameters have the following meaning. 

First line: 46.15D0 - the electron and positron beam energy, EBEAM. 

Second line: 24.D-3, 58.D-3, 0.DO - the electron minimum and maximum scattering 
angle (in radians) and the minimum visible energy (in GeV), T1MIN, T1MAX, E1MIN. 
These cuts are to be interpreted as “fiducial” cuts within which the events are generated, 
before going through the triggering routine. 

Third line: the same for the positron, T2MIN, T2MAX, E2MIN. 

Fourth line: 0.5D0, 1. D-2, O.DO, O.DO, O.DO- inputs that may be required by the 
cutting routines for the triggers. These value are stored in the vector CAL0INPUT(5) via 
the common block COMMON/CALOS. In this particular input the first figure will represent the 
calorimetric threshold and the second one the half-opening angle of the electromagnetic 
jet, while the other three will go unused. 

Fifth line: 1 - flag for symmetric cuts, ISIM. The user has to specify if the experimental 
cuts asked for are ( 1 ) or not ( 0 ) symmetric for electron-positron exchange. If they are, 
choosing 1 saves computing time. 

Sixth line: 1 - flag for choosing the triggering routine, ICALO. Possible choices are: 

1. - TRIGGER1 (sample trigger) 

2. - TRIGGER2 (free slot) 

Seventh line: 1. D4, O.DO, 0, ’ SABSPV. OUT ’ - these are inputs related to the Monte 
Carlo integration and to the management of the output. Namely, the total number of 
events, EVTS, the relative accuracy limit aimed at (the program stops when the statistical 
error reaches this level), ACCLIM. the restarting flag, IRESTART (if 1 the program tries to 
restart execution from the indicated output file, if 0 it reinitializes it), and the output file 
name, OUTFILE. 


7 Test Run Output 


The input file given in Section |] produces the following output: 


CM ENERGY 


92.30000 


THETA_EL_MIN 

THETA_EL_MAX 

THETA_P0S_MIN 

THETA_P0S_MAX 

EN_EL_MIN 

EN_P0S_MIN 

CAL0 INPUTS 


.02400 

.05800 

.02400 

.05800 

.00051 

.00051 

.50000 


.01000 


.00000 


.00000 


.00000 
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SIMM. FLAG 


1 


TRIGGER CHOICE = TRIGGER1 

CUTOFF = .0004758500 


FULL (LL) 

= 

174.34870 +/- 

.35262 

0(ALPHA) (GTGT) 

= 

-8.85671 

+/- 

.06548 

BORN PART (GTGT) 

= 

182.98221 

+/- 

.29020 

H.0. GTGT AND LL 2,8 

= 

.22320 

+/- 

.21581 

VIRTUAL + SOFT (GTGT) 

= 

-91.12947 

+/- 

.18056 

HARD PHOTON (GTGT) 

= 

263.66676 

+/- 

.65744 

ORDER ALPHA (GTGT) 

= 

172.53729 

+/- 

.70125 

COMPLETE X-SECT 

= 

172.76049 

+/- 

.76344 

COMPL. X-SECT (FACT) 

= 

172.83542 

+/- 

1.07296 


NUMBER OF EVENTS: 100000 

RANLUX INITIAL SEQUENCE: 4100 

RANLUX RESTARTING SEQUENCE: 4 1 8104045 0 

SUM VECTOR: 

.1743487004930795E+08 
-.8856709767153972E+06 
.1829822083556865E+08 
.2232019045425988E+05 
-.9112947344238646E+07 
.2636667600328065E+08 
.1727604884949618E+08 
.1725372865904168E+08 

SUM"2 VECTOR: 

.4283122616525836E+10 
.5072615600831745E+08 
.4190405492369947E+10 
.4657244067723529E+09 
.1156479410248304E+10 
.1127428205612726E+11 
.8812978799192981E+10 
.7894492407254981E+10 
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